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PREFACE 


In  a  paper  entitled  "Asymptotically  Fast  Solution  of  Toeplitz  and  Related 
Systems  of  Linear  Equations,"*  Bitmead  and  Anderson  present  an  algorithm 
for  the  inversion  of  Toeplitz  and  related  matrices.  In  this  report  we  shall 
examine  their  results  and  present  some  extensions. 


This  report,  as  well  as  Reference  1,  has  as  a  fundamental  concept  that 
of  displacement  rank  of  a  matrix.  (The  proof  of  many  of  the  results  presented 
in  this  report  and  further  results  on  displacement  rank  are  presented  in 
References  1  and  2.)  The  (+)  displacement  rank  of  a  matrix  may  be  defined 
as  follows: 


For  M  an  n  x  n  matrix  given  by 


The  (+)  displacement  rank  of  M  is  equal  to  rank  (M-M+).  (The  (-)  displacement 
rank  is  obtained  in  an  analogous  way,  except  that  M  is  obtained  from  M  by 
taking  "one  step  up  the  diagonal.")  If  M  is  a  Toeplitz  matrix,  M  is  constant 
along  diagonals  and  is,  therefore,  of  displacement  rank  less  than  or  equal  to 
two.  The  displacement  rank  may  be  used  as  a  measure  of  the  "distance"  of  a 
matrix  from  being  Toeplitz. 


1.  Bitmead,  Robert  R. ,  and  Anderson,  Brian  D.O.,  "Asymptotically  Fast 
Solution  of  Toeplitz  and  Related  Systems  of  Linear  Equations," 

Linear  Algebra  and  its  Applications,  34,  pp  103-116,  1980. 

2.  Kailath,  Thomas,  Kung,  Sun-Yuan,  and  Morf,  Martin,  "Displacement  Ranks 
of  Matrices  and  Linear  Equations,"  Journal  of  Mathematical  Analysis 
and  Applications,  68,  pp  395-407,  1979. 


A  more  direct,  and  more  important,  connection  between  displacement  rank 
and  Toeplitz  matrices  is  the  following:  The  (+)  displacement  rank  of  a  matrix 

M  is  the  smallest  number  a+  so  that  M  may  be  decomposed  as  M  L.jU.,  where 

i=l 

L.  and  U-  are  lower  and  upper  triangular  Toeplitz  matrices,  respectively. 

(This  will,  in  fact,  be  our  definition  and  we  shall  show  that  the  more  natural 
definition  given  earlier  is  equivalent.)  This  type  of  decomposition  is  of 
value  for  two  reasons.  First,  an  upper  or  lower  triangular  Toeplitz  matrix 
may  be  stored  by  storing  a  single  vector.  Second,  the  multiplication  of  two 
upper  or  two  lower  triangular  matrices  may  be  effected  by  a  single  convolution 
of  two  vectors.  Thus,  for  large  matrices  of  small  displacement  rank,  such  as 
large  Toeplitz  matrices,  the  matrices  may  be  stored  more  economically  in 
decomposed  form.  Further,  as  terms  of  the  form  UL  can  be  decomposed 
economically  in  the  form  +  L^,  multiplication  of  large  matrices 

with  small  displacement  ranks  may  be  performed  economically.  The  details  of 
obtaining  and  manipulating  these  decompositions  form  Section  One  of  this 
report  and  the  proofs  of  these  results  are  presented  in  Section  Two. 

The  reader  may  have  noted  that  forming  multiple  products  of  the  form 
(L-jU^. . .  (LnUn)  will  result  in  geometric  growth  in  the  number  of  terms  in  the 
final  decomposition.  This  appears  to  be  an  intrinsic  feature  of  such 
products.  This  unfortunate  fact  seems  to  be  the  major  drawback  to  these 
techniques  and,  unfortunately,  it  is  a  very  serious  one.  As  we  shall  see  in 
Section  Four,  the  algorithm  presented  in  Reference  1,  with  modifications  for 
symmetric  Toeplitz  matrices  has  computational  complexity  which  is  OCNlog^N) 
for  N  x  N  matrices.  Due  to  the  large  number  of  terms  in  the  decompositions  of 
multiple  products,  the  coefficient  for  this  expression  is  approximately  6300. 
It  is  possible  that  this  number  might  be  smaller  in  practice,  depending  on 
statistical  properties  of  the  matrices  to  be  inverted,  but  it  seems  likely 
that  the  algorithm  will  need  substantial  alteration  before  it  becomes 
practicable. 

This  report  is  organized  as  follows:  Section  One  contains  a  brief 
summary  of  the  principle  results  needed  to  obtain  and  manipulate  (♦)  and  (-) 


decompositions  for  matrices,  together  with  a  brief  description  of  the 
algorithm  presented  in  Reference  1.  In  Section  Two  the  results  of  Section 
One,  as  well  as  some  special  results  for  symmetric  matrices,  are  proved. 
Section  Three  is  a  detailed  exposition  of  the  algorithm  given  by  Bitmead  and 
Anderson,  specialized  for  use  in  inverting  symmetric,  positive  definite 
Toeplitz  matrices.  In  Section  Four  the  performance  of  this  algorithm  is 
analyzed. 

While  the  algorithm  of  Reference  1  is  of  questionable  interest,  the 
subject  of  Toeplitz  matrices  is  not,  and  the  techniques  presented  in 
References  1  and  2  yield  valuable  insights  into  this  subject.  Beyond 
providing  an  analysis  of  an  interesting  "failure,"  it  is  hoped  that  this 
report  can  be  of  value  as  a  collection  (and  clarification/correction)  of 
results  on  the  subject  of  displacement  rank. 


SECTION  ONE 


In  this  section  the  concept  of  displacement  rank  is  defined  and  various 
results  about  it  are  stated.  These  results  will  be  established  in  Section 
Two. 

A.  BASIC  DEFINITIONS  AND  PROPERTIES 
DEFINITION  1.1 

The  (+)-displacement  rank  of  a  matrix  M  is  the  smallest  integer  a+(M) 
such  that 

«+(M) 

H*  2  L,uf  <») 

i=l 

where  L-  and  are  lower  and  upper  triangular  Toeplitz  matrices, 
respectively.  The  (-)-displacement  rank  of  M  is  defined  as  above  with  a_ 
replacing  a+  and  L i  and  U.  interchanged. 

REMARK 

It  is  obvious  that  for  a  Toeplitz  matrix,  T,  cr+(T),  or_(T)  <  2.  This 
follows  from  the  equation  T=TA*I+I’Tu=I*T^+Tu*I,  where  T^  and  Tu  are  the  upper 
and  lower  "parts"  of  T,  respectively.  It  is  not,  perhaps,  as  obvious  that 
every  matrix  has  a  finite  displacement  rank.  This  is  not  difficult  to  see 
from  first  principles,  and  it  will  follow  immediately  from  Proposition  1.1 
below. 

DEFINITION  1.2 

For  a  row  vector  x=(x^,. . .  ,x  ),  let  x=(xn,. . .  .x^.  xT  will  denote  the 
column  vector 
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L(x)  denotes  the  lower  triangular  Toeplitz  matrix  with  first  column  x^.  U(x) 
denotes  the  upper  triangular  Toeplitz  matrix  with  first  row  x.  Let  Z  be  the 
matrix  L((0,1,0,. . . ,0))  and  let  V  be  the  matrix  U((0,1,0,. . . ,0)). 

LEMMA  1.1 

For  an  n  x  n  matrix  M,  we  have 

(a)  (ZMZ').j  =  M.^  for  i>l  and  j>l. 

=  0  for  i=l  or  j=l. 

(b)  (Z'MZ).j  =  M.+1  j+1  for  i<n  and  j<n. 

=  0  for  i=n  or  j=n 

REMARK 

We  see  from  Lemma  1.1  that  ZMZ1  is  the  displacement  of  M  "one  step  down" 
the  diagonal  and  that  Z'MZ  is  the  displacement  "one  step  up." 

PROPOSITION  1.1 

Let  M  be  an  n  x  n  matrix. 

(a)  a+(M)  =  rank  (M-ZMZ'). 

(b)  <y_(M)  =  rank  (M-Z'MZ). 

(c)  If  M  is  invertible,  then  a+(M)  =  a_(M  1 ) . 

N  N 

(d)  M-ZMZ'  =  £  xTy.iffM=  £  L(x.)U(y.). 

i=l  i=l 

N  N 

(e)  M-Z'MZ  =  £  xIyiiffM=  J 

i=l  i=l 
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REMARK 


It  is  obvious  from  (a)  and  (b)  of  the  above  that  the  displacement  ranks 
of  M  are  both  less  than  or  equal  to  the  dimension  of  M. 

DEFINITION  1.3 


Let  x=(xp. 

j 

where  z  -  = 

J 

i=l 


*.xn),  y=(y1,--- ,yn). 
xiyj+i-r 


x  y  is  the  row  vector  z=(zp. . .  ,zn), 


x#y  is  the  row  vector  w=(w^,. . .  .v^)  given  by  w=(Xp. . .  ,xn,0,. . .  ,0)* 
,yn.O,... ,0).  tr(x,i)  is  the  n-vector  (xx . x. ,0,...,0). 

PROPOSITION  1.2 


Let  x=(x  ,... ,xn),  y=(y1,... ,yn). 

(a)  L(x)L(y)=L(x  y). 

(b)  U(x)U(y)=U(x  y). 

Let  ct  denote  the  jth  column  of  L(x)U(y). 

J 

Let  C.  denote  the  jth  column  of  U(y)L(x). 

J 

Let  R^  denote  the  ith  column  of  L(x)U(y). 
Let  R.  denote  the  ith  column  of  U(y)L(x). 

(c)  Cj(i)  =  x#tr(y,j) (n+j-i). 

(d)  Cj(i)  =  y#tr(x,n+l-j)(n+j-i). 

(e)  R|(j)  =  tr(x,i)#y  (n+j-i). 

(f)  R^(j)  =  tr(y,n+l-i)#x  (n+j-i). 
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Proposition  2  gives  formulas  for  products  of  triangular  Toeplitz 
matrices.  Products  of  similar  forms  are  given  by  a  single  convolution.  For 
product  of  dissimilar  forms,  an  entire  row  or  column  of  the  product  is 
obtained  by  one  convolution. 

PROPOSITION  1.3 

Let  x,  y,  Cj  and  R*  be  as  above. 

(a)  L(x)U(y)=IL(z)+U(w)I-U(x)L(y) , 

where  x  =  (0,  xR . x2),  y  =  (0,  yn,...,y2),  z  =  R*,  w  =  tr(C*.  n-1). 

(b)  U(y)L(x)  =  L(z ' ) - I  +  I-U(w')  -  L(£)U(y), 

where  x  and  y  are  as  in  (a),  z‘  =  C^,  and  w1  =  (0,R1(2),. . . .R-^n)). 

REMARK 

The  contents  of  this  result  is  that  LU  (or  UL)  products  may  be  turned 
into  UL  (respectively,  LU)  products  with  two  convolutions  and  some  rearrange¬ 
ment  of  vector  entries. 

PROPOSITION  1.4 

Let  P  be  an  n  x  n  matrix  of  rank  m,  and  let  A  be  a  nonsingular  m  x  m 
minor.  Let  xl  be  the  column  of  P  corresponding  to  the  ith  column  of  A,  let  w 
be  the  row  of  P  corresponding  to  the  kth  row  of  A,  and  let  W  be  the  m  x  n 
matrix  whose  kth  row  is  w. .  If  y.  is  the  jth  row  of  A  "Sil,  then 
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REMARK 


Propositions  1.2  and  1.4  allow  us,  in  principle  av  least,  to  form  (+)  or 
(-)  decompositions  of  arbitrary  matrices  M. 

PROPOSITION  1.5 

Let  T  be  an  invertible  matrix  subdivided  into  retangular  blocks;  T^, 


T21* 

and  T22,  i 

T  = 

All  T12 

\T21  T22 

Let  S=T  1  be  similarly  subdivided.  Then,  if  T^,  is  square  and  invertible, 
(a)  Sn  =  T^  +  T^  T12(T22  -  T21  TxJ  T12)  1  T21  T^. 


(b)  S12  =  -TjJ  T12(T22  -  T21  TxJ  T12)  1. 

(c)  s21  =  -(t22  -  t21  txJ  t12)  1  t21  T^. 

(d)  S22  =  (T22  -  T21  TxJ  T12)  1. 


PROPOSITION  1.6 


Let  T  be  as  in  Proposition  1.5.  Let  each  T-.  be  square  and  of  size 

*  J 

n  x  n.  Let  T  have  a  (♦)  decomposition  of  the  form 


IKCy^ . y^)). 


(a)  S22  has  (-)  displacement  rank  less  than  or  equal  to  aQ.  Also, 
the  T. .  have  (+)-decompositions  of  the  following  forms: 

’  J 
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where  C*^(j)  is  the  jth  entry  of  the  nth  column  of  L(x^) 
U(y^)  and  R*^(j)  is  the  jth  entry  of  the  nth  row  of  the 
same  matrix. 


B.  AN  OUTLINE  OF  THE  ALGORITHM 


The  algorithm  described  in  Bitmead  and  Anderson1,  is  a  recursive  one. 
Given  a  matrix  of  size  2n,  the  matrix  is  subdivided  and  the  algorithm  is 
applied  to  matrices  of  size  n.  These  matrices  are  recombined  to  give  the 
inverse  of  the  2n  x  2n  matrix.  More  precisely,  the  algorithm  proceeds  as 
follows: 

“o 

Assume  that  T  =  ^  L(x.)  U(y.)  is  given. 

i=l 

1.  Subdivide  T  into  T^,  T^,  T ^ ,  T ^  using  Proposition  1.6. 

“o 

2.  Use  the  algorithm  to  write  T-^  as  T^  =  ^  Uj  Lj.  (That  such  a 

j=l 

decomposition  exists  if  T-^  is  invertible  follows  from  Proposition  1.5(a)  and 
Proposition  1.1(c).) 

3.  Use  the  results  of  step  2,  Propositions  1.2,  1.3  and  1.6  to  write 


P 

(T22  ‘  T21  T11  T12)  =  I  (In  fact,  generally,  the  number 

k=! 


3  2 

of  terms  is  3cr  +  6a„  +  4a  +  2. ) 
0  0  0 


-1 


4.  Using  Propositions  1.2,  1.4  and  1.1,  express  (I22  '  ^21  ^11  T12^  as 


^  0..  (That  this  is  possible  follows  from  Proposition  1.6(a)  and 

i=l 


Proposition  1. 1(c). ) 


5. 


Use  the  algorithm  and  the  result  of  step  4  to  write  S22  as 

“o 

S22  =  (T22  “  T21  T11  T12)  1  =  S  Ui  Li' 

i=l 

6.  Using  the  result  of  step  5  and  step  2,  calculate  minors  S^,  S^2,  S2^ 
of  T  where  the  notation  is  as  in  Proposition  1.5.  Express  each  of  the 

P 

minors  as  a  sum  of  terms  Sab  =  ^  U.  L*.  (The  number  of  terms  depends  on 

i=l 

■»  p 

the  subscripts  "a"  and  "b."  For  S22,  0  =  aQ.  For  S12  and  S21,  0  =  3aQ  +  3a^. 
For  Su,  »  =  9a®  ♦  lSoJ.  <ta;j  ♦  a,.) 

Using  the  results  of  step  6  and  employing  Propositions  1.1,  1.2  and  1.4, 
write  S  as 


i=l 


(This  is  possible  by  Proposition  1.1(c).) 

In  the  next  section,  we  shall  prove  the  results  quoted  in  this  section, 
as  well  as  stating  and  proving  some  special  results  for  symmetric  matrices. 
In  the  third  section  we  shall  give  a  more  complete  exposition  of  this 
algorithm  for  the  special  case  of  positive  definite  symmetric  Toeplitz 
matrices. 
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SECTION  TWO 


In  this  section  the  proofs  of  Propositions  1.1-1. 6  are  given,  as  well  as 
some  special  results  for  symmetric  matrices. 

LEMMA  1.1 

For  an  n  x  n  matrix  M,  we  have 

(a)  (ZMZ1  )i j  =  M.^  for  i>l  and  j>l. 

=  0  for  i=l  ur  j=l. 

(b)  (Z'MZ)^ j  =  Mi+1  j+1  for  i<n  and  j<n. 

=  0  for  i=n  or  j=n. 

PROOF 

(a)  Definition  1.2,  Z..  =  6,  ,  k,Z’  =  6m  . 

IK  i  i,k  mj  m,j-l 

n 

Thus,  (ZMZ1  )| j  =  £  6..1>k  M^  6m  ._v  which  implies  (a). 
k,m=l 

n 

(b)  ( Z 1 MZ ) -j j  =  £  6.  ^  M^  6m,i  y  which  implies  (b). 

k,m-l 

LEMMA  2.1 

Let  x=(x1,. . . ,xn),  y=(y1,. . . ,yn).  If  we  define  x$,  yt=0  for  s,  t<0,  then 

(a)  tUx)  UO)]^  y^... 

m=l 

(b)  [U(y)  L(x)]u  =  £Wk  V„.r 

m'=l 
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PROOF 


From  Definition  1.2, 

[L(x)]jt  =  xjn_,  and  [U(y)]kj  =  yjn_k. 

with  the  convention  that  xgl  yt=0  for  s,  t<0.  From  these  two  relations,  (a) 
and  (b)  follow  immediately. 

LEMMA  2.2 

Let  M  be  an  n  x  n  matrix.  If  there  exist  row  vectors  a^,...,a£ 
bl*.  *  *  ,bjn  such  ^at 


£ 


i=l 


then  rank  M  <  min(rank{a.} ,  rankfb.. }).  Further,  if  rank(M)  =  m,  then  there 
exist  row  vectors  a-j.-.-.a^,  b1,...,bJJ)  such  that  the  relation  given  above 
holds  for  M. 

PROOF 

As  a  linear  operator  on  column  vectors,  the  range  of  M  is  a  subspace  of 

T  T  T 

span  a  I  and  the  range  of  M  a  subspace  of  span  b.  .  Since  rank  M  =  rank 

T  1  1 

M  ,  the  inequality  given  above  holds. 

Next,  if  rank  M  =  m,  then  there  exists  a  spanning  set  of  rows  of  M: 

b-,,...,b  If  v.  is  the  jth  row  of  M,  then  there  exist  coefficients 
i  ro  j  u — 

a^,...^^  such  that 
m 

v.  =  ]T  bj  ^or 

i=l 

Let  a^  =  (a^,. . .  ,a^).  Then  the  relation  given  above  holds  for  M. 


LEMMA  2.3 

Let  x  =  (x1,...,xn),  y  =  (y1,... ,yn),  and  let  x  =  (xn,... ,x1),  y  = 
(yn,--,yl).  Then 

(a)  L(x)  U(y)  -  Z  L(x)  U(y)  Z'  =  xT  y. 

(b)  U(y)  L(x)  -  V  U(y)  L(x)  Z  =  yT  x. 

PROOF  (a) 

[L(x)  U(y)  -  Z  L(x)  U(y)  Z']^.  = 

n  n 

2  Vl-m  *M-m  -  2  V.J'j-.' 

iifI  m=l 

by  Lemmas  1.1  and  2.1.  But  this  is  just  x.  yj. 

PROOF  (b) 

As  above 

(UL  -  Z'ULZ)^  = 

n  n 

S  V+l-k  xm'+l-£  "  (1_6k,n)(1'6£,n)  Sym'-k  xm'-£* 
m'=l  m'=l 

This  equals  yn+1_k  xn+i-£»  which  is  the  k£  entry  of  yT  x. 

LEMMA  2.4 

If  A  and  B  are  n  x  n  matrices,  rank(I-AB)  =  rank(I-BA). 

PROOF 

Since  A  and  B  are  n  x  n,  it  is  enough  to  show  that  the  dimensions  of  the 
null  spaces  of  I-AB  and  I-BA  are  equal.  However,  if  x^,...,xk  are  independent 
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and  satisfy  x^  =  ABXj,  then  clearly  Bx. ,...,Bxj  are  independent  (else  ABx. 

=  x •  are  dependent)  and  (Bx.)  =  BA(Bx-). 

1  J  J 

PROPOSITION  1.1 

Let  M  be  an  n  x  n  matrix. 

(a)  a+(M)  =  rank(M  -  ZMZ'). 

(b)  ci_(M)  =  rank(M  -  Z'MZ). 

(c)  If  M  is  invertible,  then  a+(M)  =  a_(M  *). 

n  n 

(d)  M  -  ZMZ'  =  x{  y.  iff  M=  J  L(Xi)  U(y.). 

i=l  i=l 

n  n 

(e)  M  -  Z'MZ  =  £  x{  y.  iff  M  =  J  U(i.)  L(y - ) . 

i=l  i=l 

PROOF  (d) 

The  map  A  -*•  A  -  ZAZ'  is  linear.  Also,  by  induction  on  the  dimension,  one 
can  easily  see  that  A  -  ZAZ1  =  0  implies  A  =  0.  It  is  clear  from  Lemma  2.3 
that  if  M  =^L(x.)  U(y^),  that  M  -  ZMZ'  =  ^xT  y^.  Conversely,  if  M  -  ZMZ' 
is  equal  to^xl  y. ,  then  let  =T L(x.)  U(y.).  Then  (M  -  M  )  - 

M  II  O  AW  II  U 

Z(M  -  MQ)  V  =  0,  and  so  M  =  M0> 

PROOF  (e) 

This  follows  just  as  for  (d). 
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PROOF  (a) 


Say  M  =  ^  L(x.)  U(y.j).  Then,  by  (d)  and  Lemma  2.2,  rank(M  -  ZMZ')  <  a+ 
Also,  if  rank(M  -  ZMZ1)  =  0,  then,  by  (d)  and  Lemma  2.2,  0  >  a+. 

PROOF  (b) 

This  is  immediate  from  (e)  and  Lemma  2.2. 

PROOF  (c) 

If  M  is  invertible,  using  Lemma  2.4,  a+(M)  =  rank(M  -  ZMZ')  = 
rank(I  -  (ZMXZ'M-1))  =  rank(I  -  (2,M“1)(ZM)>  =  rank(M_1  -  Z'M-1Z)  =  a.(M-1). 
(This  is  Theorem  1  of  Kailath,  Kung,  and  Morf,  [2].) 

PROPOSITION  1.2 

Let  x  =  (x1#... ,xn),  y  =  y2, —  »yn)- 

(a)  L(x)  L(y)  =  L(x  y). 

(b)  U(x)  U(y)  =  U(x  y).  Let  ct  denote  the  jth  column  of  L(x)  U(y) 

J 

Let  C.  denote  the  jth  column  of  U(y)  L(x).  Let  r|  denote  the 
ith  row  of  L(x)  U(y).  Let  R^  denote  the  ith  row  of  U(y)  L(x). 

(c)  Cj(i )  =  x  #  tr(y , j )  (n+j-i). 

(d)  C j ( i )  =  y  #  tr(x,n+l-j)  (n+j-i). 

(e)  Rt ( j )  =  tr(x,i)  #  y( n+j-i). 

(f)  R^(j)  =  tr(y,n+l-i)  #  x(n+j-i). 
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PROOF  (a) 


From  Definition  1.2,  [L(x)]^j  =  Xj+j-j-  Thus, 
n  n 

[L(x)  Uy)]i j  =  £  x.+1_k  yk+1.j  =  £  xk'  y(i+l-j)+l-k'  * 
k=l  k‘=l 

by  substituting  k'=i+l-k.  But  this  last  is  just  the  (i+l-j)  entry  of  (x  y), 
which  is  what  we  wished  to  show. 

PROOF  (b)  • 

As  in  (a) 

n  n 

tU(x)  U(y)],j  =  ^  xk+1_.  yj+1_k  =  J  xR,  y(j+i-i)+1.k«  . 
k=l  k'=l 

by  substituting  k'=k+l-i.  This  sum  above  is  the  (j+l-i)  entry  of  (x  y). 
PROOF  (c) 

From  Definition  1.3,  for  a=(a^,. . . ,an),  b=(b1>. . . ,bn), 

j 

tr(a,j)  #  b(k)  =  J  a.  bk+1_. 
i=l 

letting  a=y,  b=x,  k=n+j-i  and  using  the  fact  that  w#v=v#w,  we  get  that 

j 

*  *  tr(y,j)(n*j-i)  =  £  yt 

£=1 

j 

=  2 

1=1 

n 

=  2  xi*l-t'  =  CUx>  UCy)]ij  • 

£'=  1 

The  proofs  of  (d),  (e)  and  (f)  are  similar  to  that  given  for  (c). 
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PROPOSITION  1.3 

Let  x,  y,  C*  and  R*  be  as  above. 

(a)  L(x)  U(y)  =  I-L(z)  +  U(w)-I  -  U(£)  L(y), 
where 

x  =  (0,xn,... ,x2),  y  =  (0,yn>... ,y2), 
z  =  R*  and  w  =  tr(C*,n-l). 

(b)  U(y)  L(x)  =  L(z')-I  +  I-U(w')  -  L(5)  U(y), 
where 

x  and  y  are  in  (a),  z'=Cp  and  w'=(0,R~(2) , . . . ,R^(n)). 

PROOF  (a) 

For  i<n  and  j<n, 

[L(x)  U(y)  -  Z'L(x)  U(y)Z].j  = 

Cj(i)  -  cj+i(i+1)  =  *  *  tr(y,j)  (u+j-i)  - 

x  #  tr(y,j+l)  (n+j+l-i+1) 

=  -x  #  (0 . 0,yj+r  0 . 0)  (n+j-i), 

where  the  yj+1  occurs  in  the  j+1  place.  This  last  quantity  is  just  =-x.j+^ 
yj+1,  by  Definitions  1.2  and  1.3.  (Recall  that  xk  =  xn+1_k. )  From  this,  the 
ij  element  of  LU-Z'LUZ,  for  i<n  and  j<n,  is  given  by  "(x2,x3,. . . ,xn,0)T 
(y2,. . . ,yn,0).  The  entries  of  the  nth  row  are  just  given  by  (0,0,. . . ,1)^R*. 
The  entries  of  the  nth  column  (excluding  the  nn  element,  to  avoid  duplication) 
are  given  by  tr(C*,n)T  (0,...,0,1).  Thus  LU-Z'LUZ  is  given  by  (0,...,0,1)T  R* 
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+  tr(C*,n)T  (0,...,0,1)  -  (x2,x3,...,xn.O)T  (y2,. . . ,yn,0).  The  result  then 


follows  from  Proposition  1.1(e). 


PROOF  (b) 


This  follows  just  as  in  (a). 


PROPOSITION  1.4 


Let  P  be  an  n  x  n  matrix  of  rank  m,  and  let  A  be  a  non-singular  m  x  m 
minor.  Let  xl  be  the  column  of  P  corresponding  to  the  ith  column  of  A,  let  wk 
be  the  row  of  P  corresponding  to  the  kth  row  of  A,  and  let  W  be  the  m  x  n 
matrix  whose  kth  row  is  w^.  If  is  the  jth  row  of  A_1W,  then 

m 

p = i  • 

i=l 


PROOF 


...  J 


Let  X  be  the  n  x  m  matrix  whose  ith  column  is  x'.  Then,  what  we  wish  to 
show  is,  if  we  let  P'-XA  that  P=P ' . 


Let  us  examine  the  matrix  X. 


x  = 
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Here,  we  have  boxed  the  rows  in  X  which  correspond  to  rows  of  A.  (That  is, 
the  rows  in  P  which  contain  the  rows  of  A  are  labeled  Similarly, 

the  matrix  W  has  the  form 


where  the  m  x  m  matrix  formed  by  the  boxed  columns  is  also  A. 

It  is  immediate  from  its  definition  that  the  n  x  n  matrix  P1  has  rank<m. 
Further,  it  is  clear  from  the  above  discussion  of  X  and  W  that  the  matrices  P 
and  P’  agree  on  the  ij  elements  where  i=r,,...,rm  or  j=Cp . . . ,Cm.  We  wish  now 
to  show  that  P=P‘. 


By  interchanging  rows  and  columns,  we  may  take  P  and  P'  so  that  r.=i  and 
C.=j.  Thus,  we  have 

J 


where  P  and  P‘  are  both  n  x  n  and  of  rank  m.  Since  the  columns  of  A  are 
independent,  a  column  of  W-^  may  be  expressed  in  only  one  way  in  terms  of  the 
columns  of  A.  But  this  implies  that  the  columns  of  P^  and  P^  are  equal,  or 
P=P\ 

LEMMA  2.5 

Let  P  be  an  n  x  n  matrix  of  rank  m.  If  the  ranks  of  the  rows  indexed  by 
i  1 , '  *  *  *  i»n  ancl  of  the  co^umns  indexed  by  are  m,  then  the  minor 

(P.  •  )  k,£=l,...,m  is  invertible. 

Vi 
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PROOF 


This  result  is  clearly  independent  of  any  row  and  column  rearrangements, 
so  that  we  may  assume  that  ik=k  and  j£=jL  Let  us  denote  the  principle  m  x  m 
minor  of  our  (possibly  rearranged)  matrix  P  by  M,  and  let  us  assume  that  M  has 
rank  k.  Finally,  let  us  take  the  columns  of  P  to  have  been  so  arranged  that 
the  first  k  columns  of  M  are  a  basis  for  the  columns  of  M.  Then  P  is  of  the 
form 


where  ft  has  rank  k,  (M|R)  has  rank  m,  has  rank  m  and  P  has  rank  m.  Let  us 
suppose  that  k<m.  Then,  since  the  rank  of  (M|R)  is  m,  there  must  exist  a 
column  of  R  which  is  independent  of  the  columns  of  M.  But  then,  the 
corresponding  column  of  in  P  must  be  independent  of  (^)in  P,  and  thus  the 
rank  of  P  must  be  greater  than  m.  This  is  a  contradiction  and  so  k=m  and  M  is 
invertible. 

PROPOSITION  2.1 

Let  P  be  a  symmetric  n  x  n  matrix  of  rank  m  and  let  M  be  an  invertible 
m  x  m  minor.  Then  there  exists  an  invertible  m  x  m  minor  A,  such  that  A  is 
symmetric  with  respect  to  the  diagonal  of  P  and  the  rows  of  P  which  contain  A 
also  contain  M. 

PROOF 

The  rows  of  P  which  contain  M  have  rank  m.  This  is  also  true  of  those 
columns  of  P  obtained  by  transposing  these  rows,  since  P  is  symmetric.  The 
result  follows  from  Lemma  2.5  if  we  take  A  to  be  the  minor  formed  by  the 
"intersections"  of  these  rows  and  columns. 

REMARK 

This  result  will  be  of  use  to  us  in  modifying  the  algorithm  of  Section  1 
for  the  special  case  of  T  symmetric. 


LEMMA  2.6 


REMARK 

Lemma  2.6  is  needed  for  the  proof  of  Proposition  1.5. 
PROPOSITION  1.5 


For  T,  T. j  as  in  Lemma  2.7,  let  S  =  T  1  and  S  be  subd‘< . .  Jed  in  a  manner 
similar  to  that  of  T.  Then 


■  iliawuijBl 


PROOF 


This  follows  easily  from  Lemma  2.7  by  matrix  multiplication. 


PROPOSITION  1.6 


Let  T  be  as  in  Proposition  1.5  with  each  of  the  T.j  of  size  n  x  n. 
have  a  (+)-decomposition  of  the  form 

T  =  f  Lax^,...^5))  UKy*1*,...^)). 

i=l 


Then 


(a)  S22  has  (-^displacement  rank  <aQ.  Also,  the  T^ .  have 

(+)-decompositions  of  the  following  forms: 
a_ 


(b) 


(c) 


TU=  \  LUxJ0,...^15))  U((y^',...,y^O). 


,(i)  „(D, 


i=] 

a. 


T12  J  L((xJ1\...,xJ1^))  U((yJjJ,...,y^))  + 


i=l 

\ 


(x<f) . x<'>)  (0,y^ . y'1’)  I  -I. 


(d) 


r21  =  £  L((xnll . x2n))}  + 


1=1 


I-U 


(0,x<<> . x<’>)  (y‘<> . y<'>) 


1=1 


<•>  T22  ’  I  L«’£l . U«^l . £>»  ♦ 

1=1 


L( V) •  +  I-U(W). 


Let  T 
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Here 


v  =  ^  (C*(i)(n),...,C*(i)(2n-l))  and 

or, 

w=  y  (0,R*(i)(n+l),...,R+(i)(2n-l)), 

1=1 

where  C*^(j)  is  the  jth  entry  of  the  nth  column  of  L(x^)  U(y^)  and 
Rp^(j)  is  the  jth  entry  of  the  nth  row  of  the  same  matrix. 

PROOF  (a) 

$22  is  the  lower  right  hand  minor  of  size  n  x  n  in  the  2n  x  2n  matrix  S. 
Since  S  has  (-)-displacement  rank  aQ,  by  Proposition  1.1(c),  rank  (S-Z1 SZ)=ao. 
But  S22~7-'$22^  15  3  minor  that  matri*  and,  therefore,  (Here 

Z,Z'  are  assumed  to  be  of  the  same  size  as  the  matrices  $,$22  in  expressions 
where  they  occur. ) 

PROOF  (b) 

To  prove  (b)-(e),  it  is  clearly  enough  to  show  these  results  for  aQ=l. 
Consider  the  matrix  T=L(x1>. . . .X2n)  U(y^,. . .  »y2n)-  The  minor  is  clearly 
the  product  of  the  first  n  rows  of  L  with  the  first  n  columns  of  U,  which  is 
simply  the  expression  in  (b). 

PROOF  (c) 

Tj2  is  the  product  of  the  first  n  rows  of  L  with  the  last  n  columns  of  U. 
However,  since  the  first  n  rows  of  l  are  zero  from  the  (n+1)  entry  on,  = 
4l*U12*  where  the  notation  has  the  obvious  meaning.  Since  U12  1-5  of  the  form 
U((yn+i,".,y2n)  +  U(0,yn,. .  .yg)),  and  since  the  product  of  lower  triangular 
Toeplitz  matrices  is  lower  triangular  Toeplitz,  we  have  (c). 
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SECTION  THREE 


In  this  section  a  version  of  the  algorithm  outlined  in  Section  One  will 
be  presented  in  more  detail.  We  shall  be  concerned  with  the  special  case  of  T 
a  symmetric  positive  definite  Toeplitz  matrix.  In  describing  the  algorithm, 
we  shall  make  use  of  three  "procedures,"  which  we  shall  refer  to  as  SUBDIVIDE, 
DECOMP  and  PRODUCT.  We  give  descriptions  of  these  procedures  below,  together 
with  discussions  of  their  storage  requirements  and  computational  complexity. 
Following  the  descriptions  of  the  procedures,  we  present  the  algorithm. 

A.  THE  PROCEDURES 

(1)  SUBDIVIDE 

This  procedure  takes  a  (+)-decomposition  (of  a  2n  x  2n  matrix  T)  of  the 
form  Ux^)  U^)  +  L^)  U^)  and  returns  (+)-decomposition  of  the  4  n  x  n 
minors  described  in  Proposition  1.6.  The  lengths  of  the  decompositions  of 
Tll’  T12’  T21’  and  T22  are’  respectively,  2,  3,  3  and  4.  If  0(n)  is  the 
number  of  computations  needed  to  convolve  two  n- vectors  (here  we  make  no 
distinction  between  x*y  and  x#y),  this  procedure  requires  40(n)  computations, 
plus  some  overhead  which  we  shall  ignore. 

(Here  we  are  not  distinguishing  between  the  various  types  of  arithmetic 
operations  in  our  analysis.) 

The  space  required  to  store  the  "matrix"  T  goes  up  from  8n  to  24n  as  a 
result  of  this  procedure,  assuming  that  the  original  decomposition  of  T  is  not 
saved. 

(2)  DECOMP 

Let  R  denote  a  symmetric  positive  definite  n  x  n  matrix  with  (+)(or(-))- 
displacement  rank  <2.  This  procedure  takes  a  (+)(or(-))-decomposition  of 
arbitrary  length  N,  and  returns  a  minimal  (+)(or(-))-decomposition.  Since  the 
(+)  and  (-)  versions  of  the  procedure  are  completely  analogous,  we  shall 
examine  only  the  (+)  version. 
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Let  P=R-ZRZ‘.  By  our  hypothesis  on  R,  P  is  symmetric,  has  rank  <2,  and 
the  element  P^  is  positive,  which  implies  that  the  first  row  and  first  column 
of  P  are  both  non-zero.  Thus,  if  we  wish  to  find  a  maximal  minor  of  P,  by 
Proposition  2.1  this  minor  is  either  of  the  form  (P^)  (and  the  rank  of  P  is 
1),  or  of  the  form 


Therefore,  all  that  is  necessary  is  to  calculate  the  determinants  of  2  x  2 
minors  of  the  form  given  above.  Thus,  we  need  to  calculate  the  first  row  and 
the  diagonal  elements  of  P. 

If 

N  N 

R  =  J  L(xi )  U(yi ) ,  then  P  =  ^  x!y . , 
i=l  i=l 

and  calculating  the  first  row  of  P  requires  N*n  multiplications  and  N*n 
additions.  Similarly,  calculating  the  diagonal  of  P  takes  2N*n  computations. 
Calculation  of  the  appropriate  determinants  requires  3n  computations  and 
comparison  of  the  results  to  select  the  "best"  minor  requires  approximately 
and  additional  n  calculations.  Finally,  calculation  of  the  minimal  decompo¬ 
sition  takes  additional  6n  computations  if  P  has  rank  2  and  n  computations  if 
P  has  rank  1.  Thus,  the  total  number  of  computations  is  at  most  (2N  +  10) -n. 
The  required  storage  goes  from  2Nn  to  at  most  4n.  (In  the  event  that  this 
algorithm  or  one  like  it  is  ever  implemented,  this  procedure  will  need  to  be 
examined  for  numerical  stability.  It  seems  likely  that  this  procedure  is  the 
critical  area  of  the  algorithm  in  this  regard.) 

(3)  PRODUCT 

Let  dT,  Dj  denote  decompositions  of  (+),  respectively  (-),  type.  This 
procedure  will  have  3  versions:  (1)  3-plus,  (2)  3-minus,  and  (3)  5-minus. 

The  3-plus  version  takes  decompositions  and  D^,  of  ranks  3  and  2, 
respectively,  and  returns  a  (+)  decomposition  of  the  product  (D^)(D^)(D^) . 
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The  3-minus  version  takes  decompositions  D^  and  Dp  of  ranks  3  and  2, 
respectively,  and  returns  a  (-)-decomposition  of  the  product  (Dj)(D^)(Dj). 
Finally,  the  5-minus  version  takes  decompositions  Dp  Dp  and  D^,  of  ranks  3, 
2,  and  2,  respectively,  and  returns  a  (- ^decomposition  of  the  product 
(D^KO^HD^XDiXDj).  As  the  three  versions  are  analogous,  we  shall  examine 
only  the  3-plus  version  and  state  the  results  needed  for  the  other  two. 

In  evaluating  the  product  (D^)(D^)(D^),  there  will  be  3  x  2  x  3  =  18 
terms  of  the  form  (LjUjXl^l^XL^Up  to  decompose.  Using  Proposition  1.2  and 
2  x  18  =  36  convolutions,  these  terms  become  18  terms  of  the  form  LXUXOUj. 
Then,  by  2  x  18  =  36  more  convolutions  we  obtain,  from  Proposition  1.3,  18 
terms  of  the  form  L^(Lg*I  +  I*Ug  +  LgUg)Up  This  requires  an  additional 
4  x  18  convolutions  to  yield,  finally,  3  x  18  =  54  terms  of  the  form  L7Up 
Thus,  the  total  number  of  convolutions  required  is  36  +  36  +  72  =  144,  and 
the  storage  requirement  is  an  additional  108n. 

The  3-minus  version  requires  96  convolutions  and  space  for  an  additional 
72n  numbers. 

The  5-minus  version  requires  2016  convolutions  and  space  for  an 
additional  648  terms  of  the  form  UL,  or  1296n  additional  numbers. 

B.  THE  ALGORITHM 

1/  If 

The  algorithm  for  a  symmetric,  positive  definite,  2x2  Toeplitz  matrix 
T  is  as  follows: 

(0)  Assume  that  T  has  been  decomposed  as  +  L2U2. 

(1)  If  the  input  matrix  is  of  sufficiently  small  size,  invert  it. 
Otherwise,  call  SUBDIVIDE. 

(2)  Use  the  algorithm  to  invert  Tp  and  return  Tp  as  UL  +  U'L'. 

(3)  Use  the  results  of  SUBDIVIDE,  step  (2),  and  the  3-plus  version 
of  PRODUCT  to  form  a  (+)-decomposition  of  Tg^TpTp.  Adjoin 
the  decomposition  of  T22  given  by  SUBDIVIDE  to  form  a  (+)- 
decomposition  of  the  matrix  B^g'TgjTpTp.  (Since  B  is  a 
symmetrically  located  minor  of  the  symmetric,  positive  definite 
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(4) 

(5) 

(6) 


(7) 


matrix  T_1,  B  is  symmetric,  positive  definite  matrix.  Also,  by 
Propositions  1.6(a)  and  1.1(c),  a+(B)<2.). 

Use  DECOMP  on  the  result  of  step  (3)  to  obtain  a  minimal 
decomposition  for  B. 

Use  the  algorithm  to  find  $22=B  *. 

Use  the  3-minus  version  of  PRODUCT  to  get  a  decomposition  of 

•  Use  the  5-minus  version  of  PRODUCT  to  obtain  a 

c.  1  cc  cl  11  _ .i 

decomposition  of  TijTj2S22^21^11*  and  aPPend  TjJ  t0  9et  a 
decomposition  of  S^. 

Use  (an  obvious  modification  of)  DECOMP  on  the  results  of  steps 


(5)  and  (6)  to  find  a  minimal  (-)-decomposition  of  S=T 


-T-l 
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SECTION  FOUR 


In  this  section  we  present  an  analysis  of  the  algorithm  given  in  the 
previous  section.  For  simplicity  we  treat  all  arithmetic  operations 
identically  and  assume  that  the  displacement  ranks  of  the  minors  encountered 
in  the  execution  of  the  algorithm  are  maximal.  Thus,  the  values  given  here 
for  computational  complexity  and  storage  requirements  are  upper  bounds.  It 
seems  likely,  however,  that  these  upper  bounds  in  practice  would  be  fairly 
sharp.  Since  this  algorithm  has  not  been,  and  in  its  present  form  is  likely 
never  to  be,  implemented,  this  is  only  speculation.  In  any  event,  the 
analysis  should  demonstrate  that  this  algorithm's  demands  in  terms  of  storage 
and  computational  complexity  are  too  great  for  it  to  be  practical  in  its 
present  form. 

A.  COMPUTATION  COUNT 


Let  C(n)  denote  the  number  of  operations  needed  to  invert  a  symmetric, 
positive  definite,  Toeplitz  matrix  T  given  in  the  form  of  a  (^-decomposition 
of  size  2.  Then,  from  Section  Three, 


C(2a)  =  4*0(2a) 

+  C(2a_1) 

+  144*e(2a_1) 

+  126 •2a'1 
+  C(2a_1) 

+  (96  +  2016) •0(2a"1) 

+  4-2a_1  +  1300*2a_1  + 
2-0(2a_1)  +  65O-0(2a"1)  + 
6*2a 


step  (1) 
step  (2) 
step  (3) 

step  (4)  (B  has  58  terms) 
step  (5) 

step  (6) 
step  (7) 

(by  an  analysis  similar  to 
that  given  for  DECOMP  in  §3) 


Thus  C(2a)«2C(2a~1)  +  28OO0(2a_1)  +  1400-23-1.  Now  a  FFT  of  a  vector  of 
length  N=2Y  requires  3N*y/2  operations.  The  convolutions  of  two  vectors  of 
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length  2a  1  require  padding  both  vectors  with  zeroes  to  make  them  of  length 
2a,  two  FFTs,  2a  multiplications  and  another  FFT.  Therefore,  8(2al)  =  ^*2a*a 
+  2a.  Therefore,  C(2a)»2C(2a‘1)  +  (12,600a  +  2,100)2a.  (It  is  likely  that 
these  coefficients  could  be  reduced  if  one  required  the  matrix  T  to  be  real.) 


Assuming  that  C(l)=l,  then 
n 

C(2n)«  ^  2n'a  (12,600a  +  2,100)2a  +  2n 
a=l 

=  2n(2,100*n  +  12,600  +  1) 

=  2n(6,300*n2  +  8,400‘n  +  1)  . 

That  is,  C(N)~6,300  N(log2N)2  +  8,400  N(log2N)  +  N. 

2 

Comparison  with  techniques  which  require  2N  operations  to  invert 
symmetric  Toeplitz  matrices  show  that  the  algorithm  given  here  is  as  fast  as 
these  when  3,150(log2N)2  +  4,200(log2N)  =  N,  or  for  N  >  221«2xl06. 

B.  STORAGE  REQUIREMENTS 

It  is  obvious  that  the  maximum  amount  of  storage  is  required  by  the 
algorithm  following  the  last  execution  of  step  (6).  The  principal  amount  of 
storage  is  required  to  hold  the  decompositions  of  S^  and  S2^.  This  amount  is 
684- n  numbers,  where  T  is  an  n  x  n  matrix. 

To  conclude,  the  algorithm  presented  in  Section  Three  for  inversion  of 
symmetric,  positive  definite,  N  x  N  Toeplitz  matrices: 

p 

(1)  Has  computational  complexity  asymptotic  to  6,300  N(log2N)  . 

(2)  Has  storage  requirements  of  at  least  684'N. 
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